A new dimension to Turing patterns 
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It is well known that simple reaction-diffusion systems can display very rich pattern formation 
behavior. Here we have studied two examples of such systems in three dimensions. First we 
• investigate the morphology and stability of a generic Turing system in three dimensions and then 

the well-known Gray-Scott model. In the latter case, we added a small number of morphogen 
sources in the system in order to study its robustness and the formation of connections between 
the sources. Our results raise the question of whether Turing patterning can produce an inductive 
signaling mechanism for neuronal growth. 



PACS numbers: 82.40.Ck, 47.54.+r, 05.45.-a 
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Fifty years ago Alan Turing proposed the now famous reaction-diffusion system involving two chemicals jlj to model 



I. INTRODUCTION 



biological pattern formation, and "morphogenesis" . Since then, it has been extensively used in studying various specific 
problems in mathematical biology j|, and currently there exists a vast literature on the subject, see e.g. Ref. || and 
references therein. In Turing systems nonoscillatory (in time) bifurcations occur provided that there is a substantial 
difference in the diffusion rates of the two substances. The appearance of a spatial pattern from the destabilization 
' of a homogeneous state in the absence of diffusion is known as diffusion-driven instability. 

The validity of the Turing models for describing developmental, ecological and chemical processes has been a matter 
of debate due to the difficulty of finding and identifying the morphogens, and it was not before 1989 that the first 
observation of a non-oscillatory steady pattern was found in a real chemical reaction Q . More recently, it was pointed 
out by Kondo and Asai Q that the stripes on the skin of the fish Pomacanthus Imperator had the characteristics of 
a Turing pattern, and this was verified later by two-dimensional numerical studies of growing domains 

The real problem with Turing models is that bifurcation analysis and numerical calculations in one and two di- 
mensions tend to produce very similar patterns, for a wide variety of models. A detailed study of patterns found 
in two dimensions revealed that a cubic term favors striped patterns, while a quadratic term produces spots, the 
latter being more robust. However, it is possible to produce more complicated shapes by coupling a Turing system 
to another Turing system Q|, to a mechanochemical model ||, or to a chemotactic mechanism The resulting 
i-^J • patterns bare striking resemblance to the ones found in the corresponding biological systems. As it is not easy to 
experimentally set up the stringent conditions for the appearance of a Turing instability, namely, spatial homogeneity 
and a very large ratio of the diffusion constants, there has been a lot of research towards investigating the effect on 
' patterns due to selection of having an inhomogeneous domain |lfjf , different boundary conditions , the presence of 
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^ . possible sources of morphogens 12 1, growing domains [p|, and curvature (T^ 

Recently there has been a suggestion that the concepts of positional information and chemical signali ng p resent in 
the Turing systems could play an important role as an inductive mechanism for neuronal connections [[15] . We are 
■ particularly interested in pursuing this idea, but any specific mechanism for Turing patterning in a real tissue has to 
happen in three dimensions. Since there have been relatively few systematic studies of three-dimensional patterns, 
we have recently done some preliminary numerical calculations in three dimensions p6[ . The results are analogous 
to the patterns in two dimensions, although it is not obvious that the pattern selection should be the same as in two 
dimensions. 

In this paper we present a series of numerical calculations carried out in a cube with periodic boundary conditions. 
We have chosen two simple generic models for this purpose: a well known version of the Gray-Scott model [ |L7| , 
and a general Turing model analyzed by Barrio et al. that allows the investigation of quadratic and cubic terms 
separately. 
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II. THE MODELS 



In this paper we are going to study these two models mainly numerically in three dimensions. In its general form 
the Turing system for modelling the evolution of the concentrations of two chemicals is given as 

d t u = D V V 2 U + f(U,V) 

dtV = D v V 2 V + g{U,V) (1) 

where U = U(x,t) and V = V(x,t) are the concentrations, and Du and Dy the respective diffusion constants. The 
reactions are modeled by the functions / and g which are typically non-linear. 

As has been pointed out earlier Q in relation to the studies of two dimensional systems, there may exist a number 
of different admissible modes having the same wave number. In three dimensions, as studied here, this degeneracy 
occurs already at small wave numbers and the question of which mode will dominate the evolution becomes even more 
difficult to answer. This is obvious when one recalls that in a finite grid the admissible modes are not continuous but 
discrete, i.e., 




where L x / y / z denote the system size in respective directions and n x / y / z the respective wave number indices. 



A. General Turing system 

As our first model, we studied a general Turing system proposed by Barrio et al. [Q. The equations of motion are 
obtained by Taylor expansion around the stationary uniform solution (U c , V c ). Keeping terms up to cubic order, the 
equations of motion read as follows 

dtu = D5\I 2 u + au(l — riv 2 ) + v(l — r-iu) 

dtv = S\7 2 v + v(/3 + ariuv) + u("f + r^v), (3) 

where u = U — U c and v = V — V c making the point (u, v) = (0, 0) a stationary solution. The constant 5 is a 
scaling factor and D is the ratio between the diffusion constants of the two chemicals. It is important to notice that 
D 1 is required for the diffusion-driven instability to occur. The parameters r\ and r-i of the non-linear interactions 
correspond to, as can be verified by symmetry arguments, stripes and spot patterns, respectively. 

The linear analysis for this system was done by Barrio et al. M. For completeness, the results that are relevant for 
this study are summarized here. The system described by Eq. (ra) has a stationary solution for (u, v) = (0, 0) and a 
second solution for 

-(a + 7) 
v = — u. 

(3+1 

By setting a = —7 it is possible to enforce (0, 0) as the only stationary solution. The dispersion relation for the 
linearized equation can be found by noticing that the spatial variation of u and v can be given in a plane wave form, 
and in order to find the eigenvalues standard Floquet analysis with u = Uoexp(At) and v — VQexp(Xt) is used. This 
leads to 

A 2 - BX + C = (4) 

where B=(a + /3- Sk 2 (l + £>)), C = {a - 5Dk 2 )(f3 - 8k 2 ) + a, and k 2 = k ■ k. 

The onse t of instability is subject to the following conditions: if a > then j3 < —a, if a < then j3 < — 1, and 
a — 2y/aD > (3D. With these conditions, the wavenumber of the most unstable mode is given by 

, 2 _ D{a - /3)-(D+ 1)t/oD 



K ~ SD{D-1) ' (5) 

This equation was used to select parameters for the simulations. It is worth noticing that the algebraic forms of B 
and C in Eq. (||) make it difficult to analyze the system further. This is why Eq. (0) is very useful. Fig. [I] shows the 
real parts of the eigenvalues plotted against the wave number k for the three modes we have used. 
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FIG. 1: The dispersion relation of Eq. Q for selected modes. The parameters were D = 0.516, a — 0.899, f3 = —0.91 and 
8 = 2 for k = 0.45 (dashed line), D = 0.122, a = 0.398, = -0.4 and 8 = 2 for k = 0.84 (dash-dot line), D = 0.516, a = 0.89, 
= —0.99 and 8 = .25 for k — 0.96 (solid line). The region above A(fc) = bounds the wave number values for unstable modes. 



B. Gray-Scott model 

The second system we study is the Gray-Scott [|18| model corresponding to two irreversible chemical reactions 

U + 2V -► 3V 

V - P, (6) 

i.e., the reaction of chemical U with two parts of V produces three parts of V. Due to the irreversible nature of the 
reactions, chemical P is an inert product. To write down the differential equations for the process, it is assumed 
that chemical U is fed in the reaction with constant rate F, the inert product is removed from the system, and that 
the second reaction is described by the rate constant K . The equations of motion for the concentrations of the two 
chemicals u = u(x,t) and v = v(x,t) in dimensionless units are as follows 

d t u = D u \J 2 u-uv 2 + F(l-u) 

dtv — D v \J 2 v + uv 2 — (F + K)v, (7) 

where it has been further assumed that chemical U is being consumed by the reaction at a rate that is proportional 
to uv 2 . The diffusion coefficients for the two chemicals are D u and D v , respectively. 

The Gray-Scott model was studied analytically, and numerically in two dimensions, by Pearson [jl?! who mapped 
the phase diagram for the system in terms of the two rate constants for the reactions. This model exhibits a very 
rich behavior ranging from time-independent steady solutions to chaotic, oscillatory, and to time-dependent phase 
turbulent behavior. Furthermore, Vastano et al. ]19[ | have shown that the system develops spatially steady patterns 
even when the two diffusion constants are equal. This behavior is particular to the one-dimensional case and it has 
not been observed in other dimensions. 
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In this study, we first performed numerical simulations of the Gray-Scott model in three dimensions. After that, we 
randomly added a small number of sources for chemical V. The motivation for adding the sources was to investigate 
the robustness of patterns and the formation of connections between the sources, thus attempting to mimic the 
development of a biological network, e.g., the formation of synaptic contacts between neurons. Currently, it is unclear 
whether Turing systems can produce such an inductive signaling mechanism for neural patterning Jl5|, ttq, p0| but the 
idea is nevertheless plausible. 



III. NUMERICAL CALCULATIONS 

First, we will present results from the numerical simulations of the general Turing system as given by Eq. (j^). 
Unless otherwise mentioned, all the simulations were performed in a cubic domain of grid size 50 x 50 x 50 using 
periodic boundary conditions. The well-known Euler algorithm was used for time integration with time step dt = 0.05 
for the general Turing system and dt = 1.0 for the Gray-Scott model. 

In addition, to test the robustness of the patterns, we performed a small number of simulations using Eq. (|) with 
an added uncorrelated Gaussian noise source rj(x,t) with the first and second moments defined as (r](x,t)) = and 
(r](x, t)r](x' ,t')) = 2eS(x — x')5(t — t'). The angular brackets denote an average and e is the intensity of the noise. 



A. General Turing system 

In their study of the two-dimensional Turing system Barrio et al. |t]] pointed out that the quadratic term in Eq. (j^) 
favors spots whereas the cubic term enhances stripes. They also noticed that the spots were more robust. 

Fig. [2] shows patterns obtained from simulations of the three-dimensional system. In the top row, the parameters 
were chosen to favor lamellar patterns, i.e., T2 — 0. The figures imply that starting from completely random initial 
conditions, it is not likely for the system to converge into a purely lamellar state. This was confirmed by a number of 
extensive simulation runs. 

The difference between Fig. ||a and Fig. ||b is that the parameters a and (3 have been selected in such a way as 
to enhance modes with wavenumbers k = 0.45 and k — 0.84, respectively. These are the most unstable modes in 
the sense of the linear approximation for parameter sets a = 0.899, (3 — —0.91, S — 2, D — 0.516, and a = 0.398, 
(3 = —0.40, 6 = 2, and D = 0.122, respectively. These particular choices of parameters were made in order to facilitate 
the comparison with the 2D results. The reason for these choices is that in the first case, in the sense of the linear 
analysis, there are very few admissible modes with positive growth rates and zero wavenumber degeneracy. As it is 
clear from Eq. (^), degeneracy becomes increasingly important at higher wavenumbers. Ideally, the choice k = 0.45 
should produce patterns with wave vector k = 2Tr(n x /L,0/L,0/L), i.e., favoring strongly the lamellar phase when 
r2 = 0. However, as seen in Fig. ||a, the system has not reached a lamellar state even after 500 000 time steps but it 
has converged to a mixed state instead. The simulations confirmed that to be the stable final state. 

Fig. |2|b displays the situation where k = 0.84 is favored. In this case, there are more closely spaced admissible 
modes leading to a competition between them. It is important to notice that the system has a finite size and thus 
there are, in the sense of the linear analysis, only certain admissible modes (Eq. 

From the topmost figures of Fig. || it is clear that the competition between the modes in the three-dimensional 
Turing system can lead to very interesting morphologies. In two dimensions, as studied by Barrio et al. [Q (Figs. 2a 
and 2b therein), the corresponding patterns display stripes with a very small (or zero) number of topological defects. 
The defects can be considered as reminiscents of the more complicated pattern selection in 3D since the 2D stripe 
patterns can be seen as cuts from a 3D lamellar system. However, as seen in Figs. and b, the 3D case displays much 
richer behavior and it is very difficult to obtain a purely lamellar pattern starting from random initial conditions. It 
should also be observed that the difference between Figs. 2a and 2b in [?J could be obtained by varying only the wave 
length (6). However, in three dimensions the qualitative difference of the two quantitatively different cases is clearer 
(Figs. |a and b). 

In the bottom row of Fig. [| the parameters r\ and have been selected to favor spots, or spherical shapes (these 
can be compared directly to Figs. 2e and 2f in Barrio et al. [Q). The spherical patterns turned out to be very robust. 
In the simulations, these morphologies developed very fast and were also stable against random Gaussian fluctuations 
that were used to test the stability of the patterns against perturbations as discussed earlier. 

To study whether it is possible to obtain pure modes, we set up the system in such a way that the initial conditions 
should favour only one selected mode. Fig. [|a shows a situation where the system was prepared in such a way that 
a layer of chemical U was set only in the midplane of the simulation box to provide favorable initial conditions for 
the lamellar structure to develop. The parameters are the same as used in Fig. Ufa. In Fig. ||b chemical U was put 
only in locations corresponding to the Ill-plane in a triangular mesh, i.e., favoring hexagonal symmetry of the spots 
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FIG. 2: Patterns obtained from simulations of the general Turing model after 500 000 iterations using random initial conditions 
on every site of the lattice: uo,vo £ (0, 1). Top row: n = 3.5 and T2 = 0; bottom row: n = 0.02 and = 0.2. Left column, 
the parameters were chosen to favor k = 0.45; in right column to favor k = 0.84. 



obtained as in Fig. 0c. In both cases chemical V was initialized uniformly over the cube as in Fig. @. The stability of 
these patterns was also tested against Gaussian fluctuations and both of them turned out to be robust, although the 
planes tended to align in a different way and curve slightly when a substantial amount of noise was used. 

To explore the effects of multiple modes on the morphology of the system, we tuned the parameters in such a way 
that mode k — 0.96 was the most unstable one. In the simulations we used parameters a = 0.89, /3 = —0.99, S = 0.25, 
D = 0.516. The non-linear parameters were r\ = 0.02 and = 0.2 corresponding to spotty patterns. Fig. |^a shows 
a stabilized configuration where the pattern shows a tubular-like structure. 

To observe competition between spherical and lamellar structures we isolated mode k = 0.45 and set both r\ and 
T2 different from zero. Using non-linear parameters r\ = 3.5 and r2 = 0.2 resulted in the competition clearly seen in 
Fig. |Jb. It took 2 000 000 time steps to stabilize the structure. 



B. Gray-Scott model 

Finally, we performed simulations in two and three dimensions using the Gray-Scott model of Eq. (fr|). In these 
simulations we have used D u — 0.125 and D v = 0.05 for the two diffusion constants. By scanning the phase space we 



FIG. 3: Patterns obtained from simulation of the general Turing model after 500 000 iterations. Left: Initial conditions: 
chemical U set in midplane, b) initial conditions: chemical U placed on Ill-plane in a triangular mesh. 




FIG. 4: Patterns obtained using Eq. a) After 500 000 time steps. The isolated mode was k — 0.96 and non- linear 

parameters ri = 0.02 and = 0.2, and b) after 2 000 000 time steps, k — 0.45, n = 3.5 and T2 = 0.2. 



obtained very rich behavior. However, complexity and time dependence of some of the solutions caused difficulties 
in visualising the resulting patterns in three dimensions. The three-dimensional Gray-Scott model showed disordered 
lamellar-like and spotty-like phases. We also studied the case in which we distributed randomly four sources of 
chemical V in the system. The sources feed the chemical at a constant rate (+0.01). The shape of these sources is 
cross-like, having six branches in three dimensions to x, y and z-directions, respectively. 

As discussed before, the motivation for including the sources is to investigate if the Turing mechanism could be 
considered as a candidate for a mechanism describing neuronal growth. Should that be the case, the sources can 
be thought of as representing neurons whereas the growing dendrites must connect them. In the sense of pattern 
formation in 3D systems, one requirement is the formation of stable tubular patterns. The Gray-Scott model clearly 
produces them, see Fig. ||, and it is feasible to use it as the starting point to explore this topic further [p0| . The 
appearance of these tubular shapes seems characteristic for the Gray-Scott model whereas it was difficult to obtain 
in the case of the general Turing system. A systematic study is on the way to address these issues EQ] . 



FIG. 5: a) Pattern obtained in a 120 x 120 lattice with periodic boundary conditions using the two-dimensional Gray-Scott 
model in the presence of eight sources of morphogen V with parameters F = 0.065, K — 0.0625, D u = 0.125, D v — 0.05. The 
sources appear as "cross-like" patterns, b) Pattern obtained using the three-dimensional Gray-Scott model with four sources 
of V, F = 0.045, K = 0.065, D u = 0.125, D v = 0.05. 



IV. DISCUSSION 



In this paper, we have studied two reaction-diffusion models of Turing type in three dimensions. The first model is 
the general model obtained by Taylor expansion and the other the well-known Gray-Scott model. 

We have shown that by enhancing more than one mode we can clearly see the competition between different 
modes. This is manifested as defects in the topology of the three-dimensional structures, whereas the patterns in two 
dimensions have very few defects. In turn this is due to the difficulty of aligning planes compared with aligning lines 
in parallel. Also the spots appear to be less regular in the case where there are more modes competing. By favoring 
certain configurations with the choice of initial conditions we were able to drive the system to a purely lamellar 
structure or hexagonal symmetry. We tested the results against Gaussian noise and found the spots more robust in all 
cases. To obtain chaotic behavior and competition we used wide window of values for the wave vector and adjusted 
non-linear parameters. 

The long range goal of our three-dimensional simulations is to model real biological development, which can be 
done only in three dimensions. There are numerous examples of lamellar structures in nature, e.g. the skin and brain 
are formed out of layers. As we showed the Gray-Scott model produces complex and stable tubular structures, where 
the tubes show bifurcations. In light of the results in two dimensions we believe that the branches connect the sources 
also in three dimensions. It is known that neurons do not always form connections with the nearest neurons, but 
with neurons far away, behind other cells. We believe that Turing patterns could explain this spatial selectivity by 
providing the signaling pathways for neurotrophic factors. Simple diffusion of these substances cannot explain the 
complexity of neuronal patterning. 
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